The starting point for hydrograph analysis is to obtain the source
data. Let’s see how it goes with sample spas-zagorye.txt
discharge data for \(1956-2020\) year
range provided with grwat. This data is for Spas-Zagorye
gauge on Protva river
in Central European plane:
library(sf) # reading and manipulating spatial data
library(tidyverse) # general data wrangling
library(mapview) # interactive mapping of spatial data
library(grwat)
mapviewOptions(fgb = FALSE)
# this is path to sample data installed with grwat
path = system.file("extdata", "spas-zagorye.txt",
package = "grwat")
# for your own data just provide the full path:
# path = /this/is/the/path/to/discharge/discharge.csv
hdata_raw = read_delim(path,
col_names = c('d', 'm', 'y', 'q'),
col_types = 'iiid', delim = ' ') # read gauge data
head(hdata_raw) # see the data
#> # A tibble: 6 × 4
#> d m y q
#> <int> <int> <int> <dbl>
#> 1 1 1 1956 5.18
#> 2 2 1 1956 5.18
#> 3 3 1 1956 5.44
#> 4 4 1 1956 5.44
#> 5 5 1 1956 5.44
#> 6 6 1 1956 5.58grwat expects the date as a single column.
Therefore, in this case we have to construct it manually. For this it is
convenient to use make_date() function from
lubridate package:
hdata = hdata_raw %>%
transmute(Date = lubridate::make_date(y, m, d),
Q = q)
head(hdata)
#> # A tibble: 6 × 2
#> Date Q
#> <date> <dbl>
#> 1 1956-01-01 5.18
#> 2 1956-01-02 5.18
#> 3 1956-01-03 5.44
#> 4 1956-01-04 5.44
#> 5 1956-01-05 5.44
#> 6 1956-01-06 5.58The basic data structure for grwat separation function is a data frame (or tibble) with two columns: date (first column) and discharge (second column). The columns can have any names, not necessarily
DandQ.
Hydrograph data may contain empty values due to corrupted records or
missed observations at specific dates. A short summary of gap periods
can be obtained by gr_get_gaps() function:
gaps = gr_get_gaps(hdata)
gaps
#> # A tibble: 8 × 5
#> num start_date end_date duration type
#> <int> <date> <date> <drtn> <chr>
#> 1 1 1956-01-01 1970-04-09 5213 days data
#> 2 2 1970-04-10 1970-04-12 3 days gap
#> 3 3 1970-04-13 1979-09-09 3437 days data
#> 4 4 1979-09-10 1979-09-11 2 days gap
#> 5 5 1979-09-12 2011-05-09 11563 days data
#> 6 6 2011-05-10 2011-05-14 5 days gap
#> 7 7 2011-05-15 2020-12-26 3514 days data
#> 8 8 2020-12-27 2020-12-31 5 days gapgr_fill_gaps() function allows filling gaps in the data.
This function first expands the date sequence if some dates are missing
between the minimum and the maximum date in the table, and then trims
the sequence of dates if missing obervations are at the beginning or end
of the sequence. After expansion all missing observations are filled by
linear interpolation between nearest values.
You can limit the maximum gap extent by threshold autocorrelation
value (autocorr parameter) or expliсit number of
observations (nobserv parameter):
fhdata = gr_fill_gaps(hdata, autocorr = 0.7)
#> grwat: filled 10 observations using 5 days window for each variable
fhdata = gr_fill_gaps(hdata, nobserv = 10)
#> grwat: filled 10 observations using 10 days window for each variableHow the time window is computed in the first case? When the
autocorr parameter is used, the ACF (autocorrelation
function) is computed first, and then its values are used to obtain the
time shift (in days) during which the autocorrelation is higher or equal
to the specified value. For Spas-Zagorye data the time lag for \(ACF \geq 0.7\) is 6 days, which can be
inferred by gr_plot_acf() function:
gr_plot_acf(hdata, 0.5)A sole discharge data is enough to separate the hydrograph into quickflow and baseflow, but is not sufficient to predict the genesis of quickflow cases. Was it due to rain or thaw? To answer such questions you also need precipitation and temperature data. Ideally, these must be measured at the gauge. But often such data is not available. In this case you need to mine this data from external sources.
One of the ways to obtain the temperature and precipitation data is to use reanalyses such as ERA5. Reanalysis data is arranged as regular grids with specific resolution. In particular, the ERA5 data has \(31\) km or \(0.28125\) degrees resolution. To use such data you must tolerate the fact that none of the reanalysis grid nodes will coincide with your gauge. Instead, you have to use the data, which is either
The last two options are the most adequate. Let’s see how it can be done.
First, we need to read the basin spatial data:
# this is path to sample basin geopackage installed with grwat
path = system.file("extdata", "spas-zagorye.gpkg", package = "grwat")
# for your own data just provide the full path:
# path = /this/is/the/path/to/discharge/basin.shp
basin = st_read(path, layer = 'basin') # read basin region
#> Reading layer `basin' from data source
#> `/Users/tsamsonov/GitHub/grwat/inst/extdata/spas-zagorye.gpkg'
#> using driver `GPKG'
#> Simple feature collection with 1 feature and 7 fields
#> Geometry type: MULTIPOLYGON
#> Dimension: XY
#> Bounding box: xmin: 35.41204 ymin: 54.88195 xmax: 36.84138 ymax: 55.57005
#> Geodetic CRS: WGS 84
gauge = st_read(path, layer = 'gauge') # read gauge point
#> Reading layer `gauge' from data source
#> `/Users/tsamsonov/GitHub/grwat/inst/extdata/spas-zagorye.gpkg'
#> using driver `GPKG'
#> Simple feature collection with 1 feature and 0 fields
#> Geometry type: POINT
#> Dimension: XY
#> Bounding box: xmin: 36.62272 ymin: 55.03669 xmax: 36.62272 ymax: 55.03669
#> Geodetic CRS: WGS 84
mapview(basin) + mapview(gauge)Next, we can buffer the data on the specified distance to catch more
reanalysis data. For this we use gr_buffer_geo function,
which approximates geographic buffer of the specified radius:
basin_buffer = gr_buffer_geo(basin, 25000)
mapview(basin_buffer, col.regions = 'red') +
mapview(basin)Alternatively you can just buffer the gauge point, though it is less meaningful since you will grab reanalysis data that falls out of gauge’s basin:
gauge_buffer = gr_buffer_geo(gauge, 50000)
mapview(gauge_buffer, col.regions = 'red') +
mapview(gauge)grwat is packaged with daily reanalysis which covers the East European territory of Russia. It has a spatial resolution of \(0.75^{\circ} \times 0.75^{\circ}\), and temporal resolution of \(1\) (one) day. Data sources include:
The coverage of the reanalysis is shown below:
The reanalysis consists of two data files, each file is about 850 Mb in size:
Download this data using the link for using with grwat.
rean = gr_read_rean(
'/Volumes/Data/Spatial/Reanalysis/grwat/pre_1880-2021.nc',
'/Volumes/Data/Spatial/Reanalysis/grwat/temp_1880-2021.nc'
) # read reanalysis data
fhdata_rean = gr_join_rean(fhdata, rean, basin_buffer) # join reanalysis data to hydrological series
#> Joining, by = "Date"
head(fhdata_rean)
#> # A tibble: 6 × 4
#> Date Q Temp Prec
#> <date> <dbl> <dbl> <dbl>
#> 1 1956-01-01 5.18 -6.46 0.453
#> 2 1956-01-02 5.18 -11.4 0.825
#> 3 1956-01-03 5.44 -10.7 0.26
#> 4 1956-01-04 5.44 -8.05 0.397
#> 5 1956-01-05 5.44 -11.7 0.102
#> 6 1956-01-06 5.58 -20.1 0.032After reanalysis data are joined you can easily plot a map of the derived spatial configuration:
# plot spatial configuration
m = mapview(basin_buffer, col.regions = 'red') +
mapview(basin) +
mapview(rean$pts[basin_buffer, ], col.regions = 'black') +
mapview(rean$pts, cex = 1)
box = st_bbox(basin_buffer)
center = st_coordinates(st_centroid(basin_buffer))
#> Warning in st_centroid.sf(basin_buffer): st_centroid assumes attributes are
#> constant over geometries of x
m@map %>% leaflet::setView(center[1], center[2], zoom = 7)To be described
The basic separation procedure is provided by
get_baseflow() function. One of the most commonly used
approaches is the method by Lyne-Hollick (1979):
resbase = fhdata %>%
mutate(Qbase = gr_baseflow(Q, method = 'lynehollick'))
# quick look at the table
head(resbase, 10)
#> # A tibble: 10 × 3
#> Date Q Qbase
#> <date> <dbl> <dbl>
#> 1 1956-01-01 5.18 3.70
#> 2 1956-01-02 5.18 3.79
#> 3 1956-01-03 5.44 3.88
#> 4 1956-01-04 5.44 3.96
#> 5 1956-01-05 5.44 4.04
#> 6 1956-01-06 5.58 4.11
#> 7 1956-01-07 5.58 4.19
#> 8 1956-01-08 5.36 4.25
#> 9 1956-01-09 5.36 4.32
#> 10 1956-01-10 5.31 4.38
resbase %>%
filter(lubridate::year(Date) == 2020) %>%
ggplot() +
geom_area(aes(Date, Q), fill = 'steelblue',
color = 'black') +
geom_area(aes(Date, Qbase), fill = 'orangered',
color = 'black')
#> Warning: Removed 5 rows containing missing values (position_stack).
#> Removed 5 rows containing missing values (position_stack).The separation is mainly parameterized by smoothing parameter which
is a = 0.925 by default, and number of passes, which are
passes = 3 by default. Changing them affects the shape of a
baseflow component:
resbase = fhdata %>%
mutate(Qbase = gr_baseflow(Q, method = 'lynehollick',
a = 0.8, passes = 5))
resbase %>%
filter(lubridate::year(Date) == 2020) %>%
ggplot() +
geom_area(aes(Date, Q), fill = 'steelblue',
color = 'black') +
geom_area(aes(Date, Qbase), fill = 'orangered',
color = 'black')
#> Warning: Removed 5 rows containing missing values (position_stack).
#> Removed 5 rows containing missing values (position_stack).Let’s see how different separation methods act in comparison to each other:
methods = c("maxwell",
"boughton",
"jakeman",
"lynehollick",
"chapman")
plots = lapply(methods, function(m) {
resbase = fhdata %>%
mutate(Qbase = gr_baseflow(Q, method = m))
resbase %>%
filter(lubridate::year(Date) == 2020) %>%
ggplot() +
geom_area(aes(Date, Q), fill = 'steelblue', color = 'black') +
geom_area(aes(Date, Qbase), fill = 'orangered', color = 'black') +
labs(title = m)
})
patchwork::wrap_plots(plots, ncol = 2)
#> Warning: Removed 5 rows containing missing values (position_stack).
#> Removed 5 rows containing missing values (position_stack).
#> Removed 5 rows containing missing values (position_stack).
#> Removed 5 rows containing missing values (position_stack).
#> Removed 5 rows containing missing values (position_stack).
#> Removed 5 rows containing missing values (position_stack).
#> Removed 5 rows containing missing values (position_stack).
#> Removed 5 rows containing missing values (position_stack).
#> Removed 5 rows containing missing values (position_stack).
#> Removed 5 rows containing missing values (position_stack).Advanced separation aims at revealing the genesis of the quickflow.
Was it due to a rain or snowmelting? For this, a joint analysis of
discharge, temperature and precipitation time series is performed by
specialized algorithm, available from gr_separate()
function in grwat package.
The genetic separation of hydrograph is controlled by more than 30
parameters. The names and meaning of these parameters can be learned
thanks to gr_help_params() function:
View(gr_help_params())Since the number of parameters is large (31), they are organized as
list, which is then fed into gr_separate() function. First,
you get the params using the gr_get_params() function. The
only parameter of the function is reg, which indicates the
region for which the parameters must be extracted. After you got the
parameters, they can be changed by accessing the list elements:
# Расчленение
p = gr_get_params(reg = 'Midplain')
p$precdays = 5
p$ftrecdays = 85Next, you can separate the hydrograph:
sep = gr_separate(fhdata_rean, p)
#> grwat: data frame is correct
#> grwat: parameters list and types are OK
head(sep)
#> # A tibble: 6 × 11
#> Date Q Qbase Quick Qseas Qrain Qthaw Type Year Temp Prec
#> <date> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <int> <int> <dbl> <dbl>
#> 1 1956-01-01 5.18 NA NA NA NA NA NA NA -6.46 0.453
#> 2 1956-01-02 5.18 NA NA NA NA NA NA NA -11.4 0.825
#> 3 1956-01-03 5.44 NA NA NA NA NA NA NA -10.7 0.26
#> 4 1956-01-04 5.44 NA NA NA NA NA NA NA -8.05 0.397
#> 5 1956-01-05 5.44 NA NA NA NA NA NA NA -11.7 0.102
#> 6 1956-01-06 5.58 NA NA NA NA NA NA NA -20.1 0.032After the hydrograph is separated, it can be summarized in a set of variables:
vars = gr_summarize(sep)
head(vars)
#> # A tibble: 6 × 57
#> Year Year1 Year2 datestart datepolend PolProd Qy Qmax datemax Qygr
#> <int> <int> <dbl> <date> <date> <int> <dbl> <dbl> <date> <dbl>
#> 1 1956 1956 1957 1956-04-14 1956-05-05 21 18.5 467 1956-04-22 8.44
#> 2 1957 1957 1958 1957-03-28 1957-05-22 55 20.3 460 1957-04-08 9.34
#> 3 1958 1958 1959 1958-04-04 1958-05-16 42 27.4 537 1958-04-21 9.90
#> 4 1959 1959 1960 1959-03-28 1959-04-27 30 27.1 406 1959-04-16 10.6
#> 5 1960 1960 1961 1960-03-27 1960-04-26 30 29.4 406 1960-04-15 12.0
#> 6 1961 1961 1962 1961-03-11 1961-05-06 56 18.7 296 1961-04-10 10.5
#> # … with 47 more variables: Qmmsummer <dbl>, monmmsummer <date>, Qmmwin <dbl>,
#> # nommwin <date>, Q30s <dbl>, date30s1 <date>, date30s2 <date>, Q30w <dbl>,
#> # date30w1 <date>, date30w2 <date>, Q10s <dbl>, date10s1 <date>,
#> # date10s2 <date>, Q10w <dbl>, date10w1 <date>, date10w2 <date>, Q5s <dbl>,
#> # date5s1 <date>, date5s2 <date>, Q5w <dbl>, date5w1 <date>, date5w2 <date>,
#> # Wy <dbl>, Wgr <dbl>, Wpol2 <dbl>, Wpol1 <dbl>, Wpol3 <dbl>, Wpavs2 <dbl>,
#> # Wpavs1 <dbl>, Wpavthaw2 <dbl>, Wpavthaw1 <dbl>, WgrS <dbl>, WS <dbl>, …These functions from grwat package allow you to:
Graphical functions are based on ggplot2 graphics.
You can plot separations for selected years using
gr_plot_sep() function:
gr_plot_sep(sep, 1976) # plot single year
#> Warning: Removed 1 rows containing missing values (position_stack).gr_plot_sep(sep, c(2016, 2017)) # plot two years sequentially
#> Warning: Removed 1 rows containing missing values (position_stack).#> Warning: Removed 1 rows containing missing values (position_stack).
And also multiple years on the same layout:
gr_plot_sep(sep, 1976:1979, # plot four years on the same page
layout = matrix(c(1,2,3,4), ncol=2, byrow = T))
#> Warning: Removed 1 rows containing missing values (position_stack).
#> Removed 1 rows containing missing values (position_stack).
#> Removed 1 rows containing missing values (position_stack).
#> Warning: Removed 2 rows containing missing values (position_stack).To get the detailed description of available variables you can invoke
gr_help_vars():
View(gr_help_vars())Parameters can be statistically tested using
test_variables(df, ..., year = NULL, locale='EN') function.
Names of the parameters are passed comma-separated in place of
.... They are quoted, so you do not need to pass them as
character strings, just write their names:
gr_test_vars(vars, Qmax)
#> Warning: `select_()` was deprecated in dplyr 0.7.0.
#> Please use `select()` instead.
#> This warning is displayed once every 8 hours.
#> Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.
#> $ptt
#> $ptt$Qmax
#>
#> Pettitt's test for single change-point detection
#>
#> data: vl[vl_cmp]
#> U* = 410, p-value = 0.01599
#> alternative hypothesis: two.sided
#> sample estimates:
#> probable change point at time K
#> 15
#>
#>
#>
#> $mkt
#> $mkt$Qmax
#>
#> Mann-Kendall trend test
#>
#> data: vl[vl_cmp]
#> z = -3.3029, n = 59, p-value = 0.0009568
#> alternative hypothesis: true S is not equal to 0
#> sample estimates:
#> S varS tau
#> -506.0000000 23376.6666667 -0.2963403
#>
#>
#>
#> $tst
#> $tst$Qmax
#>
#> Sen's slope
#>
#> data: df.theil[[2]][fltr]
#> z = -3.3029, n = 59, p-value = 0.0009568
#> alternative hypothesis: true z is not equal to 0
#> 95 percent confidence interval:
#> -5.975610 -1.605263
#> sample estimates:
#> Sen's slope
#> -4
#>
#>
#>
#> $ts_fit
#> $ts_fit$Qmax
#>
#> Call:
#> mblm::mblm(formula = eval(frml), dataframe = df.theil[fltr, ],
#> repeated = FALSE)
#>
#> Coefficients:
#> (Intercept) Year1
#> 7609.154 -3.692
#>
#>
#>
#> $tt
#> $tt$Qmax
#>
#> Welch Two Sample t-test
#>
#> data: d1 and d2
#> t = 3.3226, df = 19.72, p-value = 0.003444
#> alternative hypothesis: true difference in means is not equal to 0
#> 95 percent confidence interval:
#> 59.93319 262.61602
#> sample estimates:
#> mean of x mean of y
#> 419.7857 258.5111
#>
#>
#>
#> $ft
#> $ft$Qmax
#>
#> F test to compare two variances
#>
#> data: d1 and d2
#> F = 1.2841, num df = 13, denom df = 44, p-value = 0.5169
#> alternative hypothesis: true ratio of variances is not equal to 1
#> 95 percent confidence interval:
#> 0.5797229 3.5489038
#> sample estimates:
#> ratio of variances
#> 1.284116
#>
#>
#>
#> $year
#> Qmax
#> 1970
#>
#> $maxval
#> $maxval$Qmax
#> [1] 780
#>
#>
#> $fixed_year
#> [1] FALSE
#>
#> $pvalues
#> N Variable Change.Year Trend
#> 1 1 Maximum annual discharge during seasonal flood wave 1970 -3.69231
#> M1 M2 MeanRatio sd1 sd2 sdRatio Mann.Kendall Pettitt
#> 1 419.7857 258.5111 -38.4 162.9446 143.7931 -11.8 0.00096 0.01599
#> Student Fisher
#> 1 0.00344 0.51694This is an example with three variables selected:
tests = gr_test_vars(vars, Qygr, date10w1, Wpol3)
tests$pvalues
#> N Variable
#> 1 1 Annual groundwater discharge ("baseflow") during water-resources year
#> 2 2 First date of 10-day window discharge during winter
#> 3 3 Seasonal flood runoff (with groundwater and rainwater)
#> Change.Year Trend M1 M2 MeanRatio sd1 sd2 sdRatio
#> 1 1979 0.09505 9.09443 13.93643 53.2 1.88389 2.38607 26.7
#> 2 1987 -0.20000 28-Jan 16-Jan -12.0 32.00000 32.00000 0.0
#> 3 1988 -0.08011 12.29313 8.81473 -28.3 5.65013 4.12216 -27.0
#> Mann.Kendall Pettitt Student Fisher
#> 1 0.00002 0.00000 0.00000 0.26180
#> 2 0.44019 0.69360 0.18386 0.91868
#> 3 0.00763 0.02211 0.00909 0.09900If you want to test all parameters, just skip variable names:
tests = gr_test_vars(vars)
tests$year # this is a change year detected for each variable
#> Wy Wpol2 Wgr Wpol1 Wpavs2
#> 1978 1988 1977 1988 1972
#> Wpavs1 Wpavthaw2 Wpavthaw1 WgrW WW
#> 1972 1983 1972 1978 1978
#> WgrS WS Qy datestart Qygr
#> 1977 1972 1978 1972 1979
#> datepolend PolProd Qmax datemax Wpol3
#> 1970 1988 1970 1970 1988
#> Qmmwin nommwin Qmmsummer monmmsummer Q30w
#> 1979 1965 1979 2000 1976
#> date30w1 Q30s date30s1 Q10w date10w1
#> 1995 1979 2000 1979 1987
#> Q10s date10s1 Q5w date5w1 Q5s
#> 1979 2000 1979 1987 1979
#> date5s1 Qmaxpavthaw datemaxpavthaw CountThaws DaysThawWin
#> 2000 1977 1995 1970 1985
#> Qmaxpavs datemaxpavs CountPavs DaysPavsSum CvWin
#> 1994 2002 1977 1968 1983
#> WinProd CvSum SumProd
#> 1968 1970 1968Long-term changes are tested against breaking year, which is
calculated for each variable using Pettitt test. However, if you want to
use a fixed year, you should pass the desired breaking year into
change_year parameter:
tests = gr_test_vars(vars, Qmax, Qygr, change_year = 1987)
tests$ft # Fisher F tests to compare two variances
#> $Qmax
#>
#> F test to compare two variances
#>
#> data: d1 and d2
#> F = 1.2841, num df = 13, denom df = 44, p-value = 0.5169
#> alternative hypothesis: true ratio of variances is not equal to 1
#> 95 percent confidence interval:
#> 0.5797229 3.5489038
#> sample estimates:
#> ratio of variances
#> 1.284116
#>
#>
#> $Qygr
#>
#> F test to compare two variances
#>
#> data: d1 and d2
#> F = 0.62336, num df = 20, denom df = 37, p-value = 0.2618
#> alternative hypothesis: true ratio of variances is not equal to 1
#> 95 percent confidence interval:
#> 0.2970791 1.4352702
#> sample estimates:
#> ratio of variances
#> 0.6233639Interannual changes are visualized using gr_plot_vars()
function. Its syntax is similar to gr_test_vars() and
gr_plot_sep():
gr_plot_vars(vars, Qmax) # plot one selected variablegr_plot_vars(vars, datestart) # plot one selected variablegr_plot_vars(vars, date10w1, Wpol3) # plot two variables sequentially
gr_plot_vars(vars, Qmax, Qygr, date10w1, Wpol3, # plot four variables in matrix layout
layout = matrix(c(1,2,3,4), nrow=2, byrow=TRUE)) You can add the results of statistical tests to the plot by
specifying tests = TRUE in the function call. In that case
the subtitle with test results are added, Theil-Sen slope and Pettitt
test breaking year are drawn as solid (\(p
\leq 0.05\)) or dashed (\(p >
0.05\)) lines:
gr_plot_vars(vars, date10w1, Wpol3, DaysThawWin, Qmaxpavs,
tests = TRUE) # add test information#> Warning: Removed 5 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 row(s) containing missing values (geom_path).
#> Warning: Removed 5 rows containing missing values (geom_rect).
Alternatively, you can pass to tests the result of
test_variables(), if you need to precompute it with
specific parameters (for example, by setting exclude and
year:
gr_plot_vars(vars, date10w1, Wpol3, DaysThawWin, Qmaxpavs,
tests = gr_test_vars(vars, date10w1, Wpol3, DaysThawWin, Qmaxpavs, exclude = 1990)) # add test informationBeware that in that case you should test the variables in the same order as used for plotting. If you plot variables A, B, C and supply tests for variables X, Y, Z, they will be added without any warnings, and it is your responsibility to keep them in correspondence with each other.
Finally, you can plot all variables by not supplying column
names to plot_variables() function. In that case tests (if
you want to plot them too) should also be calculated for all
variables:
gr_plot_vars(vars, tests = TRUE)Long-term changes are the differences between summarized statistics
of one variable calculated for two selected periods. Because these
statistics reflect the differences in distributions of parameters, grwat
visualizes them as box plots using gr_plot_periods()
function. The syntax is similar to gr_plot_vars() except
that you must provide either tests or year
parameter. If both are supplied then tests is prioritized
(you can also supply a fixed year when testing variables:
gr_plot_periods(vars, Qy, year = 1978)gr_plot_periods(vars, Qy, tests = TRUE)gr_plot_periods(vars, Qy, tests = gr_test_vars(vars, Qy, year = 1985))Multiple plots can be combined on one page using layout
parameter:
gr_plot_periods(vars, Qy, Qmax, date10w1, Wpol3,
tests = TRUE,
layout = matrix(1:4, nrow = 2))To plot long-term changes for all variables just skip variable names in function call:
gr_plot_periods(vars, tests = TRUE)There is also a small helper function that plots a histogram of minimal discharge month for summer and winter periods:
gr_plot_minmonth(vars, year = 1985)There is a dedicated gr_plot_matrix() function that plot
the whole time series in a convenient matrix form. Currently it supports
plotting the value, seasons and components of the discharge:
gr_plot_matrix(sep, type = 'value')
#> Warning: Removed 147 rows containing missing values (geom_raster).gr_plot_matrix(sep, type = 'season')
#> Warning: Removed 147 rows containing missing values (geom_raster).gr_plot_matrix(sep, type = 'component')
#> Warning: Removed 147 rows containing missing values (geom_raster).You can also limit the range of years that will be plotted:
gr_plot_matrix(sep, years = 1980:1995, type = 'component')
#> Warning: Removed 36 rows containing missing values (geom_raster).Sometimes you want to compare hydrographs for selected years. There
are multiple ways to do so in R. Here we will explore three of those:
overlayed, ridgeline and horiaon plots. Since these require a small
number of manipulations via existing R packages, there is no dedicated
functionality for this in grwat.
The simplest (but not necessarily most effective) way to compare multiple hydrographs is to just to combine them on one plot:
library(ggplot2)
library(lubridate)
#>
#> Attaching package: 'lubridate'
#> The following objects are masked from 'package:base':
#>
#> date, intersect, setdiff, union
sep_sel = sep |>
filter(Year %in% c(1989, 2012))
ggplot(sep_sel, aes(ymd(20000101) + yday(Date), Q,
fill = factor(Year), group = factor(Year))) +
geom_area(alpha = 0.8, position = "identity") +
geom_line() +
scale_x_date(date_labels = "%b", date_breaks = "1 month") +
theme_minimal() +
labs(x = 'Date', y = 'Discharge', fill = 'Year')It can be seen that even in the case of two hydrographs the plot is not easy to read. It becomes nearly impossible if we increase the number of years:
sep_sel = sep |>
filter(Year %in% c(1960, 1965, 1989, 2001, 2012))
ggplot(sep_sel, aes(ymd(20000101) + yday(Date), Q,
fill = factor(Year), group = factor(Year))) +
geom_area(alpha = 0.8, position = "identity") +
geom_line() +
scale_x_date(date_labels = "%b", date_breaks = "1 month") +
theme_minimal() +
labs(x = 'Date', y = 'Year')If we are interested more in the shape and relative dimensions of the
hydrographs, and do not want the exact Y axis, then we can proceed with
ridgeline plots available at ggridges
package. grwat contains a convenient
wrapper called gr_plot_ridge():
gr_plot_ridge(sep, years = c(1960, 1965, 1989, 2001, 2012))Another alternative to the ridgeline plot is the horizon plot, which
can be built by ggHoriPlot
package. Horizon plots are more compact and allow viewing more
hydrographs at once. grwat contains a convenient
wrapper called gr_plot_hori():
gr_plot_hori(sep, years = 1960:1980)To render HTML report just pass separation and variables to
gr_report() function, and provide the output file name:
report = paste(getwd(), 'Spas-Zagorye.html', sep = '/')
gr_report(sep, vars, output = report)
browseURL(report)See report generated by this command.